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We present lattice calculations for the ground state energies of tritium, helium-3, helium-4, 
lithium-6, and carbon- 12 nuclei. Our results were previously summarized in a letter publication. 
This paper provides full details of the calculations. We include isospin-breaking, Coulomb effects, 
and interactions up to next-to-next-to-leading order in chiral effective field theory. 
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I. INTRODUCTION 



Lattice effective field theory combines the theoretical framework of effective field theory 
with numerical lattice methods. In contrast with most other ab initio methods, systematic 
errors are all introduced at the beginning when defining the truncated low-energy effective 
theory. The errors can be clearly identified as either missing operators in the lattice action, 
finite volume effects, or errors from finite Euclidean-time extrapolation. Future studies 
can build upon existing calculations in a straightforward manner by including the missing 
operators, increasing the volume, or improving the Euclidean-time extrapolation. 



Lattice effective field theory has been used to study nuclear matter 



i] 



and neutron matter 



-|7|. The method has also been applied to nuclei with A < 4 using effective field theory 



with and without pions 
found in Ref. 



- llOj . A review of lattice effective field theory calculations can be 



111 ] . Reviews of chiral effective field theory can be found in Ref. [12Ml5] . 



In this paper we present the first lattice results for lithium-6 and carbon- 12 using chiral 

dng 



effective field theory. We also present the first lattice calculations to include isospin-brea 
and Coulomb effects. Our results were previously summarized in a letter publication 
This paper provides full details of the calculations. We begin by describing the lattice 
interactions in chiral effective field theory appearing at leading order, next-to-leading order, 
and next-to-next-to-leading order. This is followed by a discussion of isospin-breaking and 
Coulomb interactions. After this all unknown operator coefficients are fit using low-energy 
scattering data. We then compute the energy splitting between the triton and helium-3. We 
discuss the auxiliary-field Monte Carlo projection method and an approximate universality of 
contributions from higher-order interactions in systems with four or more nucleons. This is 
followed by lattice results for the ground state energy of helium-4, lithium-6, and carbon-12. 

II. LEADING ORDER 

The low-energy expansion in effective field theory counts powers of the ratio Q/A. Q 
is the momentum scale associated with the mass of the pion or external nucleon momenta, 
and A is the momentum scale at which the effective theory breaks down. At leading order 



(LO) in the Weinberg power-counting scheme 



17 



181 ]. the nucleon- nucleon effective potential 



contains two independent contact interactions and instantaneous one-pion exchange. As in 
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previous lattice studies we make use of an "improved" leading-order action. This improved 
leading-order action is treated completely non-perturbatively, while higher-order interactions 
are included as a perturbative expansion in powers of Q/A. 

In our lattice calculations we use the improved LO3 lattice action introduced in Ref. [f| 
with spatial lattice spacing a = (100 MeV) -1 = 1.97 fm and temporal lattice spacing 
at = (150 MeV)- 1 = 1.32 fm. We take the parameter values g A = 1.29, f % = 92.2 MeV, 
ma = m n o = 134.98 MeV. For the nucleon mass we use m = 938.92 MeV. Many of the 
calculations presented in this paper have never been attempted before, and our choice of 
spatial lattice spacing is made to optimize the efficiency of the Monte Carlo lattice calcula- 
tions. While 1.97 fm is much larger than lattice spacings used lattice QCD simulations, we 
should emphasize that we are not probing the quark and gluon substructure of nucleons but 
rather the distribution of nucleons within nuclei. Our lattice spacing corresponds with a 
maximum filling density of more than three times normal nuclear matter density. In future 
studies the same systems will also be analyzed using smaller lattice spacings. 

Throughout this discussion we first present the interactions in continuum notation and 
then later give the corresponding lattice operator. For the continuum notation we give 
matrix elements for incoming and outgoing two-nucleon momentum states. In the following 
q denotes the t-channel momentum transfer. We use r to represent Pauli matrices in 
isospin space and a for Pauli matrices in spin space. The interactions correspond with the 
amplitude, 



We use a Euclidean-time transfer-matrix lattice formalism. The transfer matrix is the 
normal-ordered exponential of the lattice Hamiltonian, : exp(— HAt): , where At equals 
one temporal lattice spacing. We use the lattice notation adopted in several previous 
publications and which is summarized in the appendix. Let V5 =0j j=i be the lattice density- 




(1) 



3 



density correlation for the spin-singlet isospin-triplet channel in momentum space, 



32 

1 \ v a^,a/-\ a 1 , a/ -k 1 \ v at ,a/-\ a\a / -\ 



(2) 



/ s,i 
Let Vs = i j=o be the density-density correlation for the spin-triplet isospin-singlet channel, 



Vs=ij= (q) = | : p at 'W V (-<7l : +^ : E Pl*®pf\-*) : 



7 5,/ 

We use these functions to write the leading-order transfer matrix, 



M LO =: exp ^ -H hcc a t - E /(?) [£5=0,1=1^=0,7=1 (g) + ^5=1,1=0^8=1,1=0(9)] 



+ 



8/X 



(4) 



Si,S2,I ni,ti<2 

The momentum-dependent coefficient function /(g) is given by 



/(?) = /o 1 ex P 



where 



cxp 



cos gj j 



cosgj, 



(5) 



(6) 



We use the value b = 0.6, which gives approximately the correct effective range for the two 
5"-wave channels when Cs=o,i=i and Cs=i,i=o are tuned to the physical S"-wave scattering 
lengths. 



III. NEXT-TO-LEADING ORDER 

At next-to-leading order (NLO) the two-nucleon effective potential includes seven contact 
interactions carrying two powers of momentum, corrections to the two LO contact interac- 
tions, and the leading contribution from the instantaneous two-pion exchange potential 



(TPEP) 



19 



23] 



^nlo = V LO + AV<°> + W + V™. 



(7) 
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The tree-level amplitudes for the contact interactions are 

A (AV {0) ) = AC + ACpta ■ t b 

and 

A = C.q 2 + C 2 k 2 + (C 3 q 2 + C A k 2 ) (a A ■ a B ) 



+ ^2 + °b) ■ [q x k 

+ C 6 (q- a A ) (q- B B ) + C 7 [a A ■ k ) ( B B ■ k ) . 
The amplitude for the NLO two-pion exchange potential is 



24 



25) 



(8) 



(9) 



384^ 



L(q) 



1ml (5& ~ tf A - 1) + g 2 (23,1 - W9 2 a ~ l) + ^f^ 2 



3<?1 



64tt 2 /; 



-L(g) [(<f • o^) (q- Sb) ~ q 2 (&A 4 &b 



where 



Ti \ 1 / a 2 ^ 2 i + q 2 + q 



4m 2 + q 2 

(10) 

(11) 



In the lattice calculations we use a low-cutoff modification of the usual power counting 
scheme. For nearly all q < A we can expand the NLO two-pion exchange potential in 
powers of q 2 j(Am 2 A. This expansion fails to converge only for values of q near the cutoff 
scale A m 2.3m n , where the effective theory already breaks down due to large cutoff effects. 
In Fig. fll]) we show the various functions appearing in the two-pion exchange potential and 
comparsions with their analytic expansions up to 0(q 2 ) and 0(q A ). We show the function 
L(q), the dimensionless 2m n pole function, 

D ^ - 4r?' < 12 > 

as well as the dimensionless function 2m n A(q). The function A(q) appears later in our 
discussion, Eq. ( 126|) . in connection with the NNLO two-pion exchange potential. In each 
case the analytic expansion approximates the full function quite well for q less than 200 MeV. 
For our chosen lattice spacing, this covers the entire range of validity expected for the low- 
energy effective theory. 

Instead of retaining the full non-local structure of V^ P q P at this lattice spacing, we simply 



use 



v LO = + V OPEP , 



(13) 
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FIG. 1: The functions L(q), D2Tr(q), and 2m 1T A(q) appearing in the two-pion exchange potential 
and comparsions with their analytic expansions up to Oiq 2 ) and 0(q 4 ). 

V NLO = V LO + AV^ + V^. (14) 

Terms with up to two powers of q from the momentum expansion of V^lo are absorbed as 
a redefinition of the coefficients in AV^ ) and V^ 2 \ 

At next-to-leading order the lattice transfer matrix is 



M NLO = M LO - a t : [AV + AV P + V q 2 + V P ^ + V S 2, q i 

The corrections to the leading-order contact interactions are 

AV = -AC : ^p at ' a (n)p at ' a (n) :, 



(15) 



AVj^IaCp-.^Pi^Pi'" 



n 



(16) 
(17) 



HJ 



and the seven independent contact interactions with two derivatives are 

V q2 = -lc q >:J2p af ' a (n)^P af ' a (n)-, 



ri.l 



(18) 
(19) 
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= -\c s >,e ■■ Y.P a s' a WiP a s ,a ^ : ' ( 2 °) 

n,S,l 

Vs\i\<? = ■ Yl Psf( H )tfPs,i a ( H ) : > ( 21 ) 

n,S,I,l 

V {q .sr = \c {q . s? : £ £ A 5 p^(n) £ A 5 ,p«>(n) :, (22) 

^/ 2 ,(g-s) 2 = ^C , / 2 i(g . s)2 : A 5 p^} a (n) A s >p a s;}(n) :, (23) 



^ S ), = -\cl=l s) . k r- :J2 *W [K Ka (n)^vpf a {n) + </(n)A,/-(n)] : 

+ 7 : E £ W \nf/(n)^i>P a sf(n) + < s }(n)A,pf ' a (n)l : 1 . (24) 



4 

n,l,S,l',I 



The densities, current densities, and symbols A/ and vf , are defined in the appendix. The 
^(igx^-fc term eliminates lattice artifacts in the spin-triplet even-parity channels. This is 
accomplished by projecting onto the isospin-triplet channel. 



IV. NEXT-TO-NEXT-TO-LEADING ORDER 

At next-to-next-to-leading order (NNLO) there are no additional two-nucleon contact 
interactions. The two-pion exchange potential contains a subleading contribution, 

<* ^NN P L E P ] = ~^Mq) ( 2m l + O i 2m l " C 3) - 

_ ti C ^T B ) A{q) {4ml + q2) [( _, ( ^ _ q2 ^ &b)] (25) 

where 

f q 

A(q) = — arctan . (26) 

However our low-cutoff expansion in powers q 2 / (4m^) reduces the NNLO two-pion exchange 
potential to a sum of contact interactions with at least four powers of q. So in this scheme 
there are no additional contributions to the two-nucleon potential at NNLO. The only new 
contributions at NNLO are due to three-nucleon interactions, 

Vnnlo — Vnlo + ^nnlo- (27) 



/ \ / 












Jt 


Jt 



(a) 



(b) 



(c) 



FIG. 2: Three-nucleon forces at NNLO. Diagrams (a), (b), and (c) show the contact potential, 



(37V) 



contact ' 



, one-pion exchange potential Vopp j an d two-pion exchange potential Vr 



(37V) 



TPE • 



Few-nucleon forces in chiral effective field theory beyond two nucleons were introduced in 



Ref. UM. In Ref. 



26| it was shown that three-body effects first appear at next-to-next-to- 



leading order (NNLO). The NNLO three-nucleon effective potential includes a pure contact 



potential, Contact > one-pion exchange potential, Vqpe > and a two-pion exchange potential, 



(37V) 



Vr 



(37V) 



TPE > 



(37V) 



V 



(37V) 



(37V) 



\NL() 1 cnhicl 1 Ol'li + ^TPE^ (28) 

The corresponding diagrams are shown in Fig. |2j 

Similar to our continuum notation for two-nucleon interactions, we write the tree-level 
amplitude for three-nucleon interactions with nucleons A, B, C. We sum over all permu- 
tations P(A,B,C) of the labels, and <ja, <?b, qc are defined as the differences between final 



and initial momenta for the respective nucleons. The amplitudes for v}^l ct and Vq PE ' are 



r (37V) 



27, 



A 



V 



(37V) 



contact 



-E 



E 



A 



V 



(37V) 



OPE 



9 A 



D 



E 



P(A,B,C) 
QA ■ °A 

q\ + ml 



{qA ■ o B ) {ta ■ r B ) 



(29) 
(30) 



Following the notation in Ref. 



28 



P(A,B,C) 

we define dimensionless parameters ce and Cd, 



and take A v = 700 MeV. 



E 



D 



f 2 A ' 



(31) 



For convenience we separately label three parts of the two-pion exchange potential, 



Vn 



(37V) 



TPE 



V 



(37V) 



TPE1 



+ V (3N) + V 



(37V) 
TPE3- 



(32) 



S 



The corresponding amplitudes are 



A 



(3N) 



TPE1 



C3 9A_) \^ 

ft \2U I Z - 



a a) (qb ■ &b) 



A 



(3N) 



TPE2 



P(A,B,C) 
2 



Iciml f 9a_\ 



{q\ + m 2 n ) (q% + m\ 

{qA ■ &a) (qb ■ &B 



(qa ■ c[b) {ta ■ T B ) 



ft WJ £ (q 2 A + K) (q B + 



(33) 



(34) 



A 



V 



(W) 



TPE3 



C4 / 9A 



E 



U ■ &a) (Qb ■ °b) 



(Qa + m D (Qb + m l\ 



Ma x q B ) ■ o c \ \{t a x t s ) ■ t 



c 



(35) 



P(A,B,C) 

The constants ci, 03,04 parameterize the coupling of the nucleon to two pions. These 
have been determined from fits to low-energy pion-nucleon scattering data, and the values 
ci = -0.81 GeV-\ c 3 = -4.7 GeV" 1 , c 4 = 3.4 GeV" 1 are used here \m,\30\. 
At next-to-next-to-leading order the lattice transfer matrix is 



M NN lo = M NLO - a t : 



T /(3iV) , V (3N) , T/ (3iV) T/ (3iV) xr(3JV)' 
•'contact 1" •'OPE ~r V TPE1 "r )/ TPE2 V TPE3 



Mi 



LO 



(36) 



From the constraints of isospin symmetry, spin symmetry, and Fermi statistics, there is 
only one independent three- nucleon contact interaction 28j, |3l| . For our lattice action the 



contact interaction V^^ct * s a product of total nucleon densities, 



•contact g -^contact ■ / 4 \H \ lh ) 



(37) 



The one-pion exchange potential Vqpe can De written as 
V§£ = -tfgS|£r E E (^(tf^Asir'^nt)) : pf^plfi^^n) : . (38) 



'OPE 



2 ^ Qn R,S,IH',S' 



The three two-pion exchange terms Vj^ v V^^ 2 , V^p^ 3 are 



2 ^,2 



T/(3JV) _ n (3A0 \^ \^ \^ 

l/ TpE1 - ^TPEl 4 j2 g 2 



n,S,In',S' n",S" 



x (A 5 «7Tj(n" ,nt)Asirj(n,n t )) : pp a j(n)pp,j(n')p a ^ a (n) : 



(39) 



9 



-2 „,2 



T/(3iV) _ n (3N) 2 9A a t 

TPE2 — -^TDT?0''^ 



'TPE2 //t 7r 



EEE 



n,7" H',S'n",S" 



(A S /7rj(n',n t )n7rj(n, n t )) 



x (A S //7rj(n // ,n t )n7rj(n,nt)) : pp}(n)p a s fj(n")p at ' a (n) : 



1/(3*0 _ n (37V) #ia t 2 



TPE3 - -^TPES a f2n 2 



n,Si,S 2 ,S 3 /i,/ 2 A n',S' n",S" 

x (A S 'n' h (n' ,n t )A Sl 'K' h {n,nt)) (A S "7r' h {n" , n t )A S2 n' h (n, n t )) 

a' ,a / ->/\ a* ,a /-»/M a^, a / -»\ 



(40) 



(41) 



The relations between these lattice operator coefficients and the coefficients in Eq. (I29H35P 

are 



n (3JV) _ _op _ 3c £ n (3AT) 
-^contact — °- C/ — X4 A ' -^OPE 



CD 



f*K 



(3JV) _ C3 
J - 



TPE1 



(37V) 
TPE2 



2d 
■72 



(3/V) _Q_ 
TPE3 " 



(42) 
(43) 



V. ISOSPIN BREAKING AND THE COULOMB INTERACTION 



In this study we include isospin-breaking terms and the Coulomb interaction. 



Isospin 

breaking (IB) in effective field theory has been addressed in the literature |32H39| . In the 
counting scheme proposed in Ref. [39|, the isospin-breaking one-pion exchange interaction 
and Coulomb potential are considered to be the same size as 0(Q 2 ) corrections at NLO. For 
the isospin-symmetric interactions we used the neutral pion mass, m n = m n o. Therefore 
the isospin-violating one-pion exchange interaction due to pion mass differences is 

A ^opep, i B] = -{^£j [( Tl ) A (n) B + (r 2 ) A (r 2 ) B ] 
x {q ■ °a) (q- ob) 



1 



q 2 + ml ± q 2 + ml _ 



(44) 



We treat the Coulomb potential in position space with the usual «em/^ repulsion between 
protons, 

A[v ^ ] = 2^(l+Jl) (l+p.) . (45) 



r 



B 



10 



However on the lattice this definition is singular for two protons on the same lattice site. 
The resolution of this problem is to include a counterterm in the form of a proton-proton 
contact interaction. For consistency we will include all possible two-nucleon contact interac- 
tions, namely, neutron-neutron, proton-proton, spin-singlet neutron-proton, and spin-triplet 
neutron-proton. Since we will fit our isospin-symmetric interaction coefficients according to 
neutron-proton scattering data, the two neutron-proton contact interactions are just linear 
combinations of the NLO interactions, AV and AVj2. This leaves two isospin-breaking 
contact interactions. In momentum space the amplitude for these contact interactions are 

l-r 3 \ fl-T 3 

A 



A(V n: 



Cn 



A(V PP ) = C ] 



pp 



2 

l + r 3 



I + T3 



(46) 
(47) 



A \ / B 

On the lattice we add these isospin-breaking terms to the NLO transfer matrix, 



M NLO — > M NLO ,ib, 



(48) 



where 



M nlo ,ib = M NLO - a t : [V OPEP ' IB + V nn + V pp ] M LO : . 



(49) 



The isospin-breaking one-pion exchange operator is 



V 



OPEP, IB 



9 A a t 



x E E E p a k%^) P i}{n 2 ) 

1=1,2 Si,S2 ni,n2 



G Sl s 2 (ni - n 2 , m n ± ) G Sl s 2 {n\ - n 2 , m^o ) 



q 7T {m w ± / 



The Coulomb interaction operator is 



V 



EM 



1 



:«em : 



E 



ni,n 2 



r(ni - n 2 ) 



\p at ' a (n 2 ) + \ P 1*{n 2 ) 



(50) 



: , (51) 



where r is the distance on the lattice. We take the value of r at the origin to be 1/2, 



r(n) = max ( -, \n\ 



(52) 



This convention choice has no observable effect since we also have a proton-proton contact 
interaction which is fitted to proton-proton scattering data. The proton-proton contact 
operator is 



V, 



pp 



o c pp- E 



1 



in) + -p^l 



71 



(53) 
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and the neutron-neutron contact operator is 



-c nn . 



1 ,t«/^ 1 ,;■ 



W - -p i=3 



( n) 



1 at ,a/->\ 



(54) 



VI. LATTICE ARTIFACTS 



In this section we discuss the relative size of lattice artifacts produced by lattice regu- 
larization. We start with lattice artifacts that break rotational invariance. Lattice reg- 
ularization reduces the full three-dimensional rotational group down to the cubic subroup. 
Lattice operators that break rotational invariance first appear at 0{Q 2 ). These include 
local two-nucleon operators with amplitude proportional to 



(55) 



and 



(r A -r B ) QiMiMr (56) 

1=1,2,3 

These operators contain terms with total spin equal to zero, two, and four. The spin-zero 
part of these operators do not break rotational invariance and are already included in our 
set of 0{Q 2 ) local operators at NLO. The spin-two and spin-four parts of these operators 
make contributions to spin-two and spin-four transition matrix elements. For example they 
generate an unphysical mixing between the 3 Si- 3 -Di channel and the 3 D 3 - 3 Gs channel. In all 
applications discussed here, however, we compute matrix elements of operators sandwiched 
between states with definite and equal values for total spin. Hence the contribution of the 
spin-two and spin-four operators must be quadratic or higher. The net result is that these 
effects appear at 0(Q 4 ). They should be included in analyses which consider corrections 
up to N 3 LO. 

In addition to local terms, there are also non-local lattice artifacts associated with the 
one-pion exchange potential. These include 0{Q 2 ) terms from the gradient coupling of the 
pion, 

(57) 



t a ■ r B 

q 2 + ml 

and the pion propagator, 



1=1,2,3 1=1,2,3 



{ta ■ t b ) (q- &a) (q- °b, 



TjI=1,2,3 it 

(q 2 + ml) 2 



(58) 
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Once again the spin-two and spin-four parts of these operators appear only quadratically 
when computing matrix elements of operators sandwiched between states with definite and 
equal values of total spin. 

The spin-zero parts of the non-local operators in Eq. (15711581) are lattice artifacts which 
break chiral symmetry. When q < m n these operators are similar to the local 0(Q 2 ) terms 
we discussed at NLO. However for m n < q < A the non-locality of these lattice artifacts 
becomes apparent. As we will see later in our discussion of 3 Si- 3 D\ mixing, there seems 
to be some signal of these artifacts in the mixing angle. The non-local 0(Q 2 ) effects can 
be removed in future lattice studies using an 0(a 2 )-improved pion lattice propagator and 
0(a 2 )-improved gradient coupling of the pion to the nucleon. Similar non-local corrections 
to the one-pion exchange potential are generated at 0(a t Q 2 /m) by the nonzero temporal 
lattice spacing. In this case, however, the effects are numerically negligible due to our small 
value for the temporal lattice spacing, a t = (150 MeV) -1 . This has been checked explicitly 
by comparing nucleon-nucleon lattice scattering data for several different temporal lattice 
spacings. 

VII. RESULTS FOR NUCLEON-NUCLEON SCATTERING 

We measure phase shifts and mixing angles using the spherical wall method (40] . This 
consists of imposing a hard spherical wall boundary on the relative separation between the 
two nucleons at some chosen radius -R wa ii- Scattering phase shifts are determined from the 
energies of the spherical standing waves, and mixing angles are extracted from projections 
onto spherical harmonics. For neutron-neutron scattering and neutron-proton scattering, 
the asymptotic radial dependence for momentum p and orbital angular momentum L is 



where R[ (r) is the radial wavefunction and Sl and Cl are Ricatti-Bessel functions of the 
first and second kind. For proton-proton scattering, however, the long-range electrostatic 
potential requires that we use Coulomb wavef unctions. We replace Si{pr) by F L (rj,pr) and 
replace C^pr) by GL{fj,pr), where 




(59) 



(60) 



V 



2p 



e- ipr c L (ri) \Fi(L + l-irj,2L + 2, 2ipr), 



(61) 
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FIG. 3: 1 Sq neutron-proton and proton-proton phase shifts versus center of mass momentum. 

(2?) 2L+1 (nr) L+1 p~ ipr T( L 4- 1 — in) 
G L (v, P r) = LJ r(2L + 2)c £ (ry) - U ( L + 1 -^L + 2,2zpr)+tF L ( V ,pr) 1 (62) 

and 

2 L e-^ 2 \T(L + l+ir))\ 
CL{T]) = f(2LT2) • (63) 

The function 1F1 is Kummer's confluent hypergeometric function of the first kind, and the 

function U is Kummer's confluent hypergeometric function of the second kind. 

In the following plots we show lattice scattering data for spatial lattice spacing a = 

(100 MeV) -1 and temporal lattice spacing a t = (150 MeV) -1 . The 1 S'o neutron-proton and 

proton-proton phase shifts are shown in Fig. [3J For comparison we show partial wave results 

from Ref . 41[ . We see that the agreement is quite good for center of mass momenta up to 



150 MeV. To constrain the neutron-neutron contact interaction, C nn , we use the neutron 



neutron scattering length, which we take to be —18 fm with an uncertainty of ±1 fm [42H45]. 
In Fig. H] we show a comparison of the 1 S'o neutron-neutron and neutron-proton phase shifts 
as calculated on the lattice. 

In Fig. [5] we plot the 3 S\ phase shift and 3 Si- 3 Di mixing angle E\ using the Stapp pa- 



rameterization 



461 ] . The agreement with the results of the Nijmegen PWA [41J for the 3 Si 



partial wave is good up to 150 MeV. The mixing angle is good at low momenta, but devi- 
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FIG. 4: Comparison of the 1 5o neutron-neutron and neutron-proton phase shifts versus center of 
mass momentum. 

ations appear at higher momenta. This discrepancy is likely due to lattice artifacts such 
as the terms previously discussed in Eq. (I57ti58j) as well as the contribution of higher-order 
interactions. In future work some improvement may be possible using an 0(a 2 )-improved 
pion lattice propagator and 0(a 2 )-improved gradient coupling of the pion to the nucleon. 
Nonetheless the physics of 3 Si- 3 Di mixing appears correct at low energies. This we can test 
by computing the quadrupole moment of the deuteron. With no additional free parameters 
to tune we find 0.22 fm 2 at leading order and 0.29 fm 2 at next-to-leading order with isospin- 
breaking contributions. The quadrupole moment is related to the strength of the mixing 
angle at low momenta. We estimate an 8% uncertainty in fitting the mixing angle in that 
regime, and so our result for the quadrupole moment with error bars is 0.29(2) fm 2 . This 
agrees well with the physical value of 0.286 fm 2 . 

In Fig. [6] we show results for neutron-proton scattering in the 1 P\, 3 Pq, 3 Pi, and 3 P2 
channels. In all cases the comparison with physical data 4l| is good up to center of mass 
momenta of 150 MeV. 



VIII. ENERGY SPLITTING BETWEEN TRITON AND HELIUM-3 

The three-nucleon system is small enough that we can use iterative sparse-matrix eigen- 
vector methods to compute energy levels on cubic periodic lattices. We fix the coefficient eg 
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FIG. 5: 3 Si neutron-proton phase shift and 3 Si- 3 -Di mixing angle versus center of mass momentum. 

as a function of Co by matching the physical triton energy at infinite volume, —8.48 MeV. 
We consider cubes with side lengths L up to 16 fm and extract the infinite volume limit 
using the asymptotic result j^ . 

E(L) = E(oo) - je~ L/Lo + O (e~^ L / Lo ^ . (64) 

The value of cd is determined from a second observable such as the spin-doublet nucleon- 
deuteron scattering phase shifts. It turns out however that the spin-doublet nucleon- 
deuteron scattering phase shift provides only a mild constraint on c^, namely that Cr> ~ 
0(1). Currently we are investigating other methods for constraining cp, including one recent 



suggestion to determine cd from the triton beta decay rate [48|. In this analysis we simply 
use the estimate cd ~ 0(1) and check the dependence of observables upon cd- 

Although the triton energy at infinite volume is used to set the unknown coefficient eg, 
the energy splitting between helium-3 and the triton is a testable prediction. The energy 
difference between helium-3 and the triton is plotted in Fig. [7] as a function of cube length. 
We show several different asymptotic fits using Eq. (16"4"|) and different subsets of data points. 
To the order at which we are working there is no dependence of the energy splitting upon 
the value of cd- Our calculations at next-to-next-to- leading order give a value of 0.780 MeV 
with an infinite-volume extrapolation error of ±0.003 MeV. To estimate other errors we take 
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FIG. 6: Pi, 3 Po) 3 Pi ; an d neutron-proton phase shifts versus center of mass momentum. 

into account an uncertainty of ±1 fm in the neutron scattering length and a 5% relative 
uncertainty in our lattice fit of the splitting between neutron-proton and proton-proton 
phase shifts at low energies. Our final result for the energy splitting with error bars is then 
0.78(5) MeV. This agrees well with the experimental value of 0.76 MeV. 



IX. HIGHER-ORDER INTERACTIONS 



In this analysis we include all operators up to next-to-next-to-leading order. Some resid- 
ual error is expected from omitted higher-order interactions starting at 0(Q 4 ). The size of 
the error depends on the momentum scale probed by the physical system of interest. For 
well-separated low-momentum nucleons no significant deviation should occur. For two nu- 
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FIG. 8: Sketch of the relative contribution of omitted operators at 0(Q A ), 0(Q 5 ), 0(Q 6 ) for 
various nucleon configurations. The relative contribution is dominated by the last case where four 
nucleons are close together. 

cleons in close proximity the systematic error should also remain very small. The properties 



of the deuteron and soft nucleon-deuteron scattering are both accurately reproduced [10]. 
For three nucleons in close proximity the error increases a bit more, and for a tight cluster 
of four nucleons it increases further. We stop at four nucleons since a localized collection 
of five or more nucleons with no relative orbital angular momentum is forbidden by Fermi 
statistics. The expected trend for systematic errors is sketched qualitatively in Fig. |HJ 

As the sketch suggests, the relative contribution is likely dominated by the last case where 
four nucleons are close together. If this hypothesis is correct then the contribution of higher- 
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order operators to low-energy phenomena should be approximately universal. Different 
higher-order operators produce roughly the same effect on low-energy data. This situation 
is analogous to the difficulty one finds in resolving the value of cp from low-energy three- 
nucleon data. One useful consequence of this universality is that most of the residual error 
can be cancelled by adjusting the coefficient of an effective four-nucleon contact term, 



This effective four-nucleon contact interaction should not be confused with the four-nucleon 
contact interaction that appears at 0(Q 6 ). We are not suggesting a rearrangement of power 
counting in chiral effective field theory. We are simply taking advantage of the expected 
universality of missing higher-order interactions. Later in our discussion we present results 
which test and appear to confirm this universality hypothesis. 

The inclusion of V^2 ive provides an opportunity to resolve another related issue that 
was noted in earlier lattice calculations. Let |4iV one _ S ite) be a configuration of four nucleons 
on a single lattice site, 

|4iV one _ si te) = al t0 (n)a\ t0 (n)al A (n)a\ A (n) 1°) • ( 66 ) 

The potential energy of this configuration is dependent upon the three-nucleon contact 
operator and the local part of the three-nucleon one-pion-exchange interaction, 

(4iV one _ S it e | V^. ontact |4iV OIie _ S it c ) = 4Z) contact , (67) 

(4iV one . site | V$g |4iV one _ sitc ) = U^D§^ ( 68 ) 

If -DcontLt or -^ope i s sufficiently large and negative, a clustering instability can be produced 
in systems with four or more nucleons. This is a lattice artifact that appears on coarse 



lattices 49j . and is similar to the clustering instability found with point-like two-nucleon 
contact interactions [Sj. That problem was solved by using improved lattice actions with 
operator smearing. An analogous technique could be adopted for the three-nucleon interac- 



tions. In Ref. lOj, however, a different approach was used. In that analysis the temporal 
lattice spacing was adjusted to ensure that the size of the cutoff-dependent three-nucleon 
operator coefficients were small. 
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In this study we use a simpler and more direct technique. Let us define 
n' {AN) -(AN \V {3N) +V (3N) + V (4N) UN • \ 

-^effective ~~ V*- 1 v onc-sitc | •'contact ' v OPE ' v effective l^ 1 v one-site/ 

= U3SSL + UD83^ E G«(3) + on,- (69) 

The problem is that the local three-nucleon terms induce an effect much the same as a four- 
nucleon contact interaction, and quite possibly a strong four-nucleon interaction. To remedy 
this we treat -Deflective as a bare counterterm that removes the dependence on D>contLt an d 
-Dope- In the following we express all lattice results in terms of the renormalized coupling 

effective ' 



X. AUXILIARY FIELDS AND PROJECTION MONTE CARLO 

For systems with more than three nucleons, sparse-matrix calculations using the lattice 
transfer matrix are not practical at large volumes. Instead we use projection Monte Carlo 
with auxiliary fields. The auxiliary-field transfer matrix for the LO3 action requires sixteen 
auxiliary fields. One auxiliary field is associated with the total nucleon density N^N, three 
fields for the spin density N^aN, three fields for the isospin density iVViV, and nine fields for 
the spin-isospin density N^BtN . Let us define M( nt \ir' I: s, s s , sj, s s ,i) as the leading-order 
auxiliary- field transfer matrix at time step n t , 

M (rit) (vr/ ? s, s s , si, s s ,i) = : exp < -H {tcc a t - ^"L ^ Ag7i^(ra, n t )p a gf(n) 



n,S,I 



+ \^(-3C s= o,i=i ~ 3C s=li/=0 ) a t ]T s(n, n t )p a '' a {n) 

n 

+ (-3C s =o,i=i + C s =i,i=o) a t ^2 Ss ^' n t)Ps' a (n) 



4 

n,S 



+ (C s =o,i=i - SC s =i,i=o) a t ^2 S K™> n t)pf' a (n) 

n,I 

+ i V ( -C s=o,j=i ~ Cs=i,/=o) S5 ^' n t)p a sf( H ) 

n,S,I 



(70) 



We can write Mlo as the normalized integral 

/ Dti' i DsDs s Ds i Dss,i e - s - ) - s -* ) M( n ')(7r^ s, s s , sj, s s ,i) 
M LO = J - j ■ - (71) 
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where S^f 1 is the piece of the instantaneous pion action at time step n t , 

St\A) = \j2<^ n ^(n,n t ) -^J^^^tKifi + lnt), (72) 

n,I n,I,l 

and si™^ is the auxiliary- field action at time step n t , 

S t ] = \Y f s(n,n t )f-\n-n')s{n\n t )+ l - £ s s (n, n*)/" 1 ^ ~ n)s s (n' , n t ) 

n,n' n,n',S 

+ 2 ^2 s /(^>^)/ _1 (^ - n^s^n 1 \n t ) + - ^ s S j(n,n t )f~ 1 (n-^)s SJ (n\n t ), 

ri,n' ,1 n,n',S,I 

(73) 

with 

f-\n -n') = ^Y. JJfj e-^~*'\ (74) 

The contributions from NLO, NNLO, isospin-breaking, and Coulomb interactions are 
treated using perturbation theory. This is done by including external sources coupled to 
densities and current densities. Let us define 

/ Dn'jDsDssDsjDssj e^' 3 ^ M^«, s, s s , s I} s s ,i,s) 

M LO (e) = J - - f , (75) 

/ Dtt'jDsDssDsiDssj e~ s ^^ ~ Ssst 

where 

=: M^K, s, s s , sj, s SJ ) exp [u^\e) + U<j?\e)\ : . (76) 
The isospin-independent couplings are 

U {nt \e) = Y^e p {n,n t )p a '> a {n) + J2e PS (n,n t )pf a (n) + J2^ SP (n,n t )A s p at ' a (n) 

n n,S n,S 

+ ^A sPsl (n,n t )A s pi' a (n) + J2 £ vfp(^ n t)^P a '' a (^) 

n,S,S' n,l 

+ J2 £ vfps( H ' n t)rfp a s ,a (n) + ^£n ; (™, n t )U.f' a (n) + ^e ni , s (n, n t )Y[f/(n), 

n,l,S rt,l n,l,S 

(77) 
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the isospin-dependent couplings are 

U { p\e) = ^2e Pl (n,n t )pf' a (n) + ^ £ Ps ,i(ri,n t )psf(n) + ^ e AsPl (n,n t )A s pf' a (n) 

n,I H,S,I n,S,I 

n,S,S',I n,l,I 

+ £ v2 ps , / (n,n t )V^5 ) } a (n) + ^2eu u (n,n t )uf/(n) + ^ e Ui SI {n, nt)n^}(n). 

n,l,S,I ft,l,I n,l,S,I 

(78) 

All of the NLO, NNLO, isospin-breaking, and Coulomb interactions are generated by func- 
tional derivatives with respect to the external source fields. 

We extract the properties of the ground state using Euclidean-time projection. Let 
j^ffree^ ^ e a g] a ^ er determinant of free-particle standing waves in a periodic cube for some 
chosen number of nucleons and quantum numbers. Let Mgu( 4 w t> e an auxiliary-field transfer 
matrix at time step n t , 



exp 



-H itee a t + A/-Csu(4)y«f ^ s(n, n t )p at ' a [ 



n) 



(79) 



We use the operator Mgu(4w( s ) ^° se ^ U P the initial state for the lattice calculation, 

\W)) = (M su(4V )^|vI/ free >, (80) 

where t' = L to a t and L to is the number of "outer" time steps. As the notation suggests, 
the operator M^ 4 ^(s) is invariant under Wigner's SU(4) symmetry [5^]. The repeated 
multiplication by Mgu*(4w( s ) ac ^ s as an approximate low-energy filter. This part of the 
Euclidean-time propagation is positive definite for any even number of nucleons invariant 



under the SU(4) symmetry 5lM53l| 



The Euclidean-time amplitude Z(t) is defined as 

Z(t) = <*(OI(M L o)^|*(0>. (81) 

where t = L ti at and L ti is the number of "inner" time steps. The transient energy at time 
t + at/2 is calculated by taking a numerical derivative of the logarithm of Z(t), 



-E ho (t+a t /2y at _ z (t + a t) 



Z(t) 



(82) 
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The ground state energy E ^o equals the asymptotic limit of the transient energy, 

£ 0) lo= Hm E h0 (t + a t /2). (83) 



t— >oo 



We calculate Euclidean-time projection amplitudes using the auxiliary-field formalism. 
For a given configuration of auxiliary and pion fields, the contribution to the amplitude Z(t) 
is proportional to the determinant of an A x A matrix of one-body amplitudes, where A is the 
number of nucleons. Integrations over auxiliary and pion field configurations are computed 
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55|. 



using hybrid Monte Carlo. Details of the method can be found in Ref. 

The perturbative contributions from NLO, NNLO, isospin-breaking, and Coulomb inter- 
actions are computed order-by-order in perturbation theory. For the first-order perturbative 
correction to the energy, it suffices to compute operator expectation values. For general 
operator O we define the Euclidean-time amplitude, 

Z (t) = W)l (M LO ) Lt ' /2 (M LO ) Ltj2 • (84) 

The expectation value of O for \^/q) is extracted by taking the large t limit of the ratio of 
Z (t) and Z(t), 

}^W- = ^o\0\^ )- (85) 

In the appendix we show precise numerical tests of the equivalence of the auxiliary-field 
Monte Carlo formalism and the original transfer matrix formalism. 



XI. RESULTS FOR HELIUM-4 



We compute the ground state energy for helium-4 in a periodic box of length 9.9 fm. For 
j^free^ we ^ a k e fog Slater determinant formed by standing waves, 

(0| a i;j (n) |V>i) oc <J ij0 5j,i, (0| a i;j (n) |^ 2 ) oc ^,0^,0, (86) 

(0| a itj (n) \tp 3 ) oc 5<,i^,i, (0| a itj (n) \ip 4 ) oc 5 itl 5 j:0 . (87) 

This produces a state with zero total momentum and the quantum numbers of the helium-4 
ground state. For each value of the Euclidean time, t, we use 2048 processors to generate 
about 5 x 10 6 hybrid Monte Carlo trajectories. Each processor runs independent trajectories, 
and averages and stochastic errors are calculated from the distribution of results from all 
processors. 
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For the numerical extrapolation in t, we use a decaying exponential for the leading-order 
energy, 

E hQ (t) ^E 0tLO + A LO e- 5E4 . (88) 

For each of the perturbative energy corrections from NLO, isospin-breaking (IB), electro- 
magnetic (EM), and NNLO interactions we use 

AE(t) « AE + AAe~ 5E - t/2 . (89) 

The unknown parameters -Eo,lo ; ^lo ; AA, and SE, are determined by least squares fitting. 
The e~ SE4 dependence in Eq. (1551) gives the contribution of low-energy excitations with 
energy gap 8E above the ground state. The e~ SE4 ^ 2 dependence in Eq. (1591) gives the 
contribution of matrix elements between the ground state and excitations at energy gap SE. 

Given the finite interval over which we measure the Euclidean-time dependence, we expect 
some exponential dependence from other energy excitations not at energy SE above the 
ground state. In order to estimate the size of the induced systematic errors, we generate an 
ensemble of different exponential fits which include dropping the two first two data points and 
then dropping the last two data points. This gives some estimate of the spread in energies 
of contributing higher energy states. In the following we quote total extrapolation errors 
which include the uncertainty due to the stochastic errors and the effect of the distribution 
in SE. In future studies we hope to improve this process further by considering different 
initial states in order to triangulate a common extrapolated value at infinite t. 

In Fig. [9] we show the energy versus Euclidean time projection for the helium-4 ground 
state with LO, NLO, IB, EM, and NNLO interactions. The plot on the left shows the 
leading-order results and the extrapolated t oo values for the higher-order contributions 
added cumulatively. These cumulative results are shown with error bars on the right edge 
of the plot. The plot on the right shows the higher-order corrections separately. For each 
case we show the best fit as well as the one standard-deviation bound. We estimate this 
bound by generating an ensemble of fits determined with added random Gaussian noise 
proportional to the error bars of each data point and also varying the number of fitted data 



points. These results are similar to those found in Ref. 



lOj using the L0 2 action. For 



c D = 1 we get -30.5(4) MeV at LO, -30.6(4) MeV at NLO, -29.2(4) MeV at NLO with 

)(4A0 
effective 



IB and EM corrections, and —30.1(5) MeV at NNLO. When the bare interaction ' 



is held fixed, the helium-4 energy decreases 0.4(1) MeV for each unit increase in cp. 
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FIG. 9: Ground state energy for helium-4 as a function of Euclidean time projection. See text for 
details. 

Apart from direct comparisons with experimental data, an independent estimate of sys- 
tematic errors due to truncation of higher-order terms can be made by comparing the dif- 
ferences among the lattice results at each order, LO, NLO, and NNLO. One caveat here is 
that sometimes the differences can be unusually small, either by chance or due to underly- 
ing physics. For example there is only a very small difference between the LO and NLO 
energies for helium-4. This can be explained by the fact that the interactions for helium-4 
are predominantly in the S'-channels, and the improved LO3 action is already quite accurate 
for 5*-wave scattering. For helium-4 we estimate a residual error of size about 1 MeV for 
the omitted interactions. This appears consistent with the 1.8 MeV deviation between the 
NNLO result and the physical binding energy for helium-4. 

For nuclei beyond A = 4, we will test the universality hypothesis for higher-order inter- 
actions by tuning the effective four-nucleon contact interaction Infective to gi ye the physical 
helium-4 energy of —28.3 MeV. The contribution of the effective four-nucleon contact in- 
teraction to the helium-4 energy is shown in Fig. [9j 
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XII. RESULTS FOR LITHIUM-6 



We compute the ground state energy for lithium-6 in a periodic box of length 9.9 fm. 
For |\[/ frcc ^ we choose standing waves, 

(0| aij(n) |V>i) oc £1,1*7,1, (0| a itj {n) \ijj 2 ) oc 8 itl 8 jfi . (90) 

(0| a itj {n) IVs) oc cos (0| a itj {n) |^ 4 ) cx 5 ifi 8 j)0 cos (91) 

(0| a hj (n) ^5) oc 5 ifi 5 3 ,i sin (0| a hJ (n) |^ 6 ) oc 5 ij0 5 j)0 sm (92) 

This combination produces a state with zero total momentum and the quantum numbers of 
the lithium-6 ground state. For each value of t a total of about 5 x 10 6 hybrid Monte Carlo 
trajectories are generated by 2048 processors. 

In Fig. [10] we show the energy versus Euclidean time projection for lithium-6. For the 
numerical extrapolation in t we use the same decaying exponential functions in Eq. (I88H89I) . 
We show the best fit as well as the one standard-deviation bound. For cd = 1 we get 
-32.6(9) MeV at LO, -34.6(9) MeV at NLO, -32.4(9) MeV at NLO with IB and EM 
corrections, and —34.5(9) MeV at NNLO. Our error estimate due to truncation at NNLO 
is about 2 MeV. Adding the contribution of the effective four-nucleon interaction ineffective 
to the NNLO result gives —32.9(9) MeV. This lies within error bars of the physical value 
—32.0 MeV. However we expect some overbinding due to the finite periodic volume. The 
finite volume analysis in Ref. 1(J found a finite volume dependence of less than 1 MeV 
for the helium-4 ground state in a periodic box of length 9.9 fm. However a larger effect 
is expected for lithium-6 due to the larger spatial distribution of the two P-shell nucleons. 
Further calculations at varying volumes will be needed to determine this volume dependence. 

Compared with helium-4, there is a much larger difference between the LO and NLO en- 
ergies for lithium-6. This may indicate additional binding coming from the NLO corrections 
in P-wave channels. The dependence of the energy on Cr> can be analyzed in several differ- 
ent ways. When the bare interaction ^effective is held fixed, the lithium-6 energy decreases 
0.7(1) MeV for each unit increase in cjj. When the effective four-nucleon interaction is ad- 
justed according to the physical helium-4 energy, the lithium-6 energy decreases 0.35(5) MeV 
per unit increase in cd- 
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FIG. 10: Ground state energy for lithium-6 as a function of Euclidean time projection. See text 
for details. 

XIII. RESULTS FOR CARBON-12 

We compute the ground state energy of carbon-12 in a periodic box of length 13.8 fm. 
For |v|> free ^ we take the Slater determinant formed by standing waves, 

(0| a itj {n) |V>4fc+i) oc 5 ifi 5 jA f k (n), (0| a itj (n) l^fc+2) oc 5 ifi S jfi f k (n), (93) 

(0| aij(n) \ip ik+3 ) oc Si^xfkin), (0| a itj {n) |V>4fc+4) oc 5 iA 5 jfi f k (n), (94) 

where 

/o(n) = l, / 1 (n) = cos^, f 2 (n) =sin^. (95) 

This combination produces a state with zero total momentum and the quantum numbers 
of the carbon-12 ground state. For each value of t a total of 2 x 10 6 hybrid Monte Carlo 
trajectories are generated by 2048 processors. 

Fig, rm shows the energy versus Euclidean time projection for carbon-12. For Co = 1 we 
get -109(2) MeV at LO, -115(2) MeV at NLO, -108(2) MeV at NLO with IB and EM 
corrections, and —106(2) MeV at NNLO. Our error estimate due to truncation at NNLO 
is about 5 MeV. The small 2 MeV difference between NLO and NNLO results is due to a 
cancellation of several larger contributions. Adding the contribution of the effective four- 
nucleon interaction ^effective ^° ^ ne NNLO result gives —99(2) MeV. This is an overbinding 
of 7% compared to the physical value, —92.2 MeV. While this agreement as a final result 
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FIG. 11: Ground state energy for carbon-12 as a function of Euclidean time projection. See text 
for details. 

would not be bad, an overbinding of 7% is actually a reasonable estimate of the finite volume 
correction for carbon-12 in a periodic box of length 13.8 fm. If so the error at infinite volume 
would in fact be much smaller than 7%. Further calculations at varying volumes will be 
needed to measure the volume dependence. 

When the bare interaction ^effective * s ne ^ fixed, the carbon-12 energy decreases 
1.7(3) MeV per unit increase in cjy. When the effective four-nucleon interaction is adjusted 
according to the physical helium-4 energy, the carbon-12 energy decreases only 0.3(1) MeV 
per unit increase in cjj. The much reduced dependence upon on cd is consistent with our 
universality hypothesis regarding systematic errors. In three-nucleon systems the value of 
cd is difficult to resolve due to similarities of the one-pion exchange three-nucleon inter- 
action and the three-nucleon contact interaction at low energies. For systems with four 
or more nucleons, the difference between these three-nucleon interactions becomes signif- 
icant. However our universality hypothesis suggests that this difference behaves like an 
effective four-nucleon contact interaction. This explains why the dependence on cd goes 
away when we include an effective four-nucleon contact interaction tuned to the physical 
helium-4 energy. 
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XIV. SUMMARY AND COMMENTS 



In this paper we have presented several new methods and results in lattice effective field 
theory. We described the first lattice results for lithium-6 and carbon- 12 using chiral effective 
field theory. This represents a significant advance in the range of problems accessible using 
lattice effective field theory. We also detailed the first lattice calculations to include isospin- 
breaking and Coulomb interactions, and computed the energy splitting between helium-3 
and the triton. The accuracy of the lattice calculations presented here are competitive with 
recent calculations obtained using other ab initio methods. Coupled cluster calculations 
without three-nucleon interactions are accurate to within 1 MeV per nucleon for medium 



mass nuclei 



56|. Constrained-path Green's function Monte Carlo calculations generally 



have an accuracy of 1% — 2% in energy for nuclei A < 12. The most recent result for carbon- 
12 is —93.2(6) MeV using AV18 and the IL7 three-nucleon force 57j. The most recent 
no-core shell model calculation for carbon-12 with the JISP16 NN interaction considers two 
different extrapolation methods to obtain values —93.9(1.1) MeV and —95.1(2.7) MeV 



We also mention some recent lattice QCD simulations in the strong coupling limit. While 
quite different from physical nuclei, the strong coupling analog of nuclei have been simulated 
for up to twelve nucleons 59 1. 

Future lattice studies should look at probing large volumes, decreasing the lattice spacing, 
and including higher-order interactions. The computational scaling with the number of 
nucleons suggests that larger nuclei are also possible. At fixed volume we find that the 
time required by one processor to generate one HMC trajectory scales with the number of 
nucleons as A 1 ' 7 for A < 16. For carbon-12 calculations the time required by one processor 
to generate one HMC trajectory scales with volume as V 1 ' 5 . For nuclei with S = and 
1 = the average sign (e' 6 *) scales as e~ 0A1A . From this scaling data we estimate that a 
simulation of oxygen- 16 would require about 1.8 TFlop-yr. 

Lattice effective field theory should prove a useful tool for few-body calculations of nuclei 
as well as many-body calculations of neutron and nuclear matter. The method is also quite 
attractive theoretically as it uses only the general principles of effective field theory. All 
systematic errors are introduced up front when defining the truncated low-energy effective 
theory. This eliminates approximation errors tied with a specific calculational tool, physical 
system, or observable. The reduction of these errors is not necessarily easy. However they 



29 



can be clearly identified as either missing operators in the lattice action, finite volume effects, 
or errors from finite Euclidean-time extrapolation. Future studies can then improve upon 
existing calculations in a straightforward manner. 
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Appendix A: Lattice notation 

The vector n represents integer- valued lattice vectors on a three-dimensional spatial lat- 
tice, and p, q, k represent integer- valued momentum lattice vectors. I — i, 2, 3 are unit 
lattice vectors in the spatial directions, a is the spatial lattice spacing, and L is the length 
of the cubic spatial lattice in each direction. The lattice time step is a t , and rit labels 
the number of time steps. We define a t as the ratio between lattice spacings, a t = a t /a. 
Throughout our lattice discussion we use dimensionless parameters and operators, which 
correspond with physical values multiplied by the appropriate power of a. Final results are 
presented in physical units with the corresponding unit stated explicitly. 

We use a and a) to denote annihilation and creation operators. We make explicit all 
spin and isospin indices, 

a 0,0 = a f,p, a o,i = o-|>> (Al) 
Qi,o = Ol,p, 0-1,1 = Ol,ti- (A2) 

The first subscript is for spin and the second subscript is for isospin. We use 77 with 
/ = 1,2,3 to represent Pauli matrices acting in isospin space and as with S — 1,2,3 to 
represent Pauli matrices acting in spin space. For the free nucleon we use the 0(a 4 )- 
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improved lattice Hamiltonian, 



H, 



free 



49 
12m 



Y Y 4* (™) 



ft i,j=0,l 



~ 4^ Y Y Y \ a h a *>j a L (") «i j - o 



n i,j=0,l 1=1,2,3 



+ ^Y Y Y \4A n ) a ^( n+2i ) +a U H ) a ^(n-2i) 



ft i,j=0,l J=l,2,3 



n tj=0,l 2=1,2,3 



(A3) 



The eight vertices of a unit cube on the lattice is used to define spatial derivatives. For 
each spatial direction I = 1, 2, 3 and any lattice function f(n), let 

^f^) = \ Y (-ir + 7(n + ^), V=v l i + v 2 2 + V &. (A4) 
We also define the double spatial derivative along direction i, 



^f{n) = f{n + i) + f{n-i)-2f{n). 
For the three-body NNLO interactions we also use the notation 



^l,J y 2,I / 3=0,l 



1. Local densities and currents 



(A5) 



(A6) 



We define the local density, 



ij'=0,l 



(A7) 



which is invariant under Wigner's SU(4) symmetry 50] . Similarly we define the local spin 
density for S = 1, 2, 3, 

P a s' a (ri)= Y alAtyMiitai'Arf)' ( A8 ) 



w'=o,i 



isospin density for I = 1, 2, 3, 



W''=0,1 



(A9) 
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and spin-isospin density for S, I = 1, 2, 3, 

Psf(ft)= E a L'(«) N jT ^(n). (A10) 

For each static density we also have an associated current density. Similar to the defini- 
tion of the lattice derivative A; in Eq. (IA4I) . we use the eight vertices of a unit cube, 

V= v x \ + z/ 2 2 + u 3 3, (All) 

for Ux, U2, V3 = 0, 1. Let u(—l) for I — 1, 2, 3 be the result of reflecting the Z th -component of 
z? about the center of the cube, 

z7(-/) = z7+ (1 - 2v l )l (A12) 

Omitting factors of i and 1/m, we can write the / th -component of the SU(4)-invariant current 
density as 

nf '"(n) = \ Yl E + *H))«u(* + V). (A13) 

^1,^2,^3=0,1 i,j=0,l 

Similarly the / th -component of spin current density is 

= i E E + *(-Q) h«i ^(n + *?), (am) 

fi,i / 2,i / 3=0,l ij,i'=0,l 

Z th -component of isospin current density is 

<' a (^) = \ E E (- 1 )' /i+1 4i(" + *H)) Hy ^ + (A15) 

"1> "2 ,1^3=0,1 i,j,j'=0,l 

and / th -component of spin-isospin current density is 

n fe(^) = i E E (-ir i+ %(n + H-l))Ws}^[ri} Jj/ a,An + ^)- (A16) 

^1,^2,^3=0,1 i,j,i',j'=0,l 
2. Instantaneous free pion action 

The lattice action for free pions with purely instantaneous propagation is 
<SW(7Tj-) = ot t {^ + 3) ^2 Trj(n, n t )7Ti(n, n t ) - a t ^ n I (n,n t )7i I (n + l,n t ), (A17) 

n,nt,I n,nt,I,l 
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where 717 is the pion field labelled with isospin index /, and m n = m n o. It is convenient to 
define a rescaled pion field, n'j, 



q w = a t (ml + 6). 



Then 



SnTT^'i) = 7; ^2 71^(71, n t ) - — ^'i(n, n t)^'i(n + i,n t ) 



n,nt,I 

In momentum space the action is 



<?7T 



n,nt,I,l 



I,k 



1 a t sr^ , 

> COS/C; 

2 (fcr V 



(A18) 
(A19) 

(A20) 



(A21) 



The instantaneous pion correlation function at spatial separation n is 

//^ \ / /7t A /-Otj <(n,n t )7r^(0,n t ) exp [-S^] 
717(71, n t j7Tj(0, rat) ) = r — 37 — - — , (no sum on I) 



J _D7Tj exp [— jS'j, 



-ik-ft 



(A22) 



where 



Ar(fc) 



1 - — E; cos h ' 



It is also useful to define the two-derivative pion correlator, Gs 1 s 2 (n), 
Gs^in) = (A 5l 7rj(7i,7it)As 2 7rj(0,7i t )^ (no sum on /) 



j o — ^ 

^1,^2,^3=0,1 u[y 2 y 3 =o,i 



+ z7 - z/, n t )vrK0, n t 



(A23) 



(A24) 



3. Pion mass differences 

We outline the modifications that result from different masses for the charged pion and 
neutral pion. Let 



q n (m n ±) = a t (ml ± + 6), q^m^) = a t {m\ + 6). 



(A25) 



The rescaled pion fields are then 

■K' 12 (n,n t ) = y/qn(m v ±)n 1) 2(n,n t ), n' 3 (n,n t ) = v^r("V>) 7r 3 (A26) 
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The momentum-space correlators for the charged and neutral pions are 

D n (k,m n ±) = 2m 1 -, (A27) 

l-^I)E.cosfc z 

^acosV (A28) 



We can now repeat the steps in Eq. ( IA24j) to define the two-derivative pion correlators 
G Sl s 2 (n,m 7T ±) and G SlS2 (n, m w o). 

Appendix B: Precision tests 

We use the three-nucleon system as a precision test of the lattice formalism and computer 
codes. The same observables are calculated using both auxiliary-field Monte Carlo and the 
exact transfer matrix without auxiliary fields. We choose a small system so that stochastic 
errors are small enough to expose disagreement at the 0.1% — 1% level. We choose the 
spatial length of the lattice to be L = 3 lattice units and set the outer time steps L to = 
and inner time steps L t . = 4. With 2048 processors we generate a total of about 10 7 hybrid 
Monte Carlo trajectories. Each processor runs completely independent trajectories, and we 
compute averages and stochastic errors by comparing the results of all processors. 

We choose |\l> frce ^) to be a Slater determinant of free-particle standing waves where 

(0| a id (n) oc 6i i0 6jfi, (0| a id (n) \ip 2 ) oc S iA 6 jfi , (0| a id (n) \ip 3 ) oc 8 ij0 S jtl . ( B1 ) 

The quantum numbers of this state correspond with helium-3 at zero momentum. At 
leading order we find an energy of —49.72(6) MeV for the Monte Carlo calculation and 
—49.7515 MeV for the exact transfer matrix. In Table [J we compare Monte Carlo results 
(MC) and exact transfer matrix calculations (Exact) for the derivative of the energy with 
respect to each NLO coefficient. Table HT1 shows the energy shifts due to the proton-proton 
contact interaction and the Coulomb interaction, and Table II III shows the derivative of 
the energy with respect to each NNLO coefficient. The numbers in parentheses are the 
estimated stochastic errors. In all cases the agreement between Monte Carlo results and 
exact transfer calculations is consistent with estimated stochastic errors. 
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TABLE I: Monte Carlo results versus exact transfer matrix calculations for the derivative of the 
energy with respect to NLO coefficients. 



NLO energy derivatives 


MC 


Exact 


d(A d E (Tc) {t)) t 104 MeV3 l 


3.722(3) 


3.72347 


{ tZ ^ V' 10 4 MeV 3 


-4.530(6) 


-4.53590 


d(AE NLO (t)) rm9 i\i„-sr5l 

d(c q2 ) \- eV J 


— z.UOO^z J 


— z.UOooo 


8(Ag NLO (t)) r 1Q 9 M y5i 


3.052(3) 


3.05148 


0(A£ NLO (t)) r 1Q 9 M y5l 


161 (3) 


16376 


JABnloW) [10 9 MeV 5] 


5.240(5) 


5.24260 


a(A^ L o 3 (t)) [1q9 MeV5] 

a ( C ( 9 -s) 2 ) 


-1.5873(9) 


-1.58896 


yg NLO (t)) [10 9 MeV 5] 


6.833(3) 


6.83234 


l^gSLoW) [i 9 MeV 5 ] 

0(^(i9XS).feJ 


0.3356(5) 


0.33702 


e>(AE NLO (t)) [1() 9 Me y5] 

a ( C / 2 ,(i9XS)-fc) 


-0.996(2) 


-0.99656 



TABLE II: Monte Carlo results versus exact transfer matrix calculations for the energy shifts due 
to the proton-proton contact interaction and the Coulomb interaction. 



IB and EM energy shifts 


MC 


Exact 


AE pp (t) [10~ 2 MeV] 


1.937(2) 


1.94128 


AE EM {t) [10- 1 MeV] 


3.712(2) 


3.71232 
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